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INTEGRATION OF CONTROL EQUATIONS AND THE PROBLEM 
OF SMALLTIME CONSTANTS 

SUMMARY 


A system of first-order, linear differential equations with constant 
coefficients is transformed into a convenient and considerably simplified system 
of differential equations, referred to as the canonical equations. Alternatively, 
a transfer function in the form of a rational function is represented by a like 
system of canonical equations. The latter equations can be readily solved 
provided the forcing function Ej n can be conveniently represented. Furthermore, 
the canonical equations which depend upon poles (eigenvalues) of large magni- 
tude may be discarded when it is evident, as is often the case, that the contri- 
butions of their solutions to the response E QU £ of the system are negligibly 
small. The elimination of such equations enables one to integrate the system 
numerically with a larger integration step-size than that which would ordinarily 
be required. 

Often it can not be easily decided which, if any, canonical equations 
can be discarded and, furthermore, E^ n can not be simply represented and may, 
in fact, not even be predictable as in cases in which it is related to E QU £ through 
a system of nonlinear differential equations. Such cases are treated by assum- 
ing E^ n to be a linear function of time over short time intervals At; that is, E. 
is expressed as 111 

E. (t) = E. (T) + ~ [e. (T + At) - E. (T)l 
in in At ]_ in' in J 

on the interval 

T < t T + At. 

Assuming the values of all significant variables to be known at time T, 
the solution E ou ^(T + At) to the canonical equations is expressed analytically 
in terms of the unknown quantity E- n ( T + At) . The analytical expression for 
E out (T + At) is then substituted into the nonlinear differential equations, 



thereby effectively eliminating the variable E out -from the problem. The re- 
sulting set of nonlinear differential equations may be integrated numerically. 

The important advantage, derived from determining E- n analytic ally and 
then eliminating it, is the very significant increase in the numerical integra- 
tion step-size frequently made possible, especially when the given system of 
linear differential equations has eigenvalues of large magnitude (or, 
equivalently, small time constants). 

As a by-product of the theory, one obtains the relationship of the 
eigenvalues and eigenvectors of a linear system of differential equations to the 
poles and coefficients of the corresponding transfer function. 

Two illustrative examples of problems with small time constants are 
given, and subroutines for implementing the theory on digital computers are 
described. 


INTRODUCTION 


This paper is concerned with the numerical integration of systems of 
equations that have the following form : 

y k -e k (t ' y >’ y *’--" y v E in , E out ) (k-l. 2,..., n d ) (la) 

E ln = *<t. yi, y* y v v E „ut> <lb) 

L(E out )/L(E ln )=p(s)/g(s) (lc) 


Differentiation is with respect to time t, L(f) signifies the Laplace transform 
F(s) of the function f(t) , and p(s) and g(s) are polynomials in s such that the 
degree of p is equal to or less than the degree of g. 

Equations (i) occur frequently in the simulation of control systems. 

The subsystem composed of (la) and (lb) represents a system of nonlinear 
differential equations such as the equations of motion of a space vehicle, 
whereas the transfer function p(s)/g(s) in (lc) represents a system of linear 
differential equations with constant coefficients which define a response function 
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E ou t<t) i n terms of a forcing function Ei n (t). (It is common to represent 
filters, actuators, and some other components of control systems by means 
of transfer functions. ) 

Very often control engineers are surprised to find that numerical integra- 
tion of equations (i) requires a prohibitively small integration step-size even 
though the dependent variables do not appear to vary rapidly enough to warrant 
such a small step. In an earlier paper [ i] the author discussed a frequent 
cause of such difficulty and showed how the problem could be treated. The 
present paper formulates the theory from a point of view more familiar to the 
engineer and more general in application. Much use has been made of experience 
gained from the extensive use of techniques previously developed [ i] . The theory 
may be applied to advantage even when a small integration step-size is not re- 
quired. Furthermore, it is applicable to the case in which the coefficients of 
the linear equations vary slowly with time. 

Other papers dealing with very similar problems are cited in the 
Reference section [2-12] . Reports [2-4] describe analytical solutions which re- 
quire elaborate coupling of first- and/or second -order linear differential 
equations. The method of Steinman [3] has the advantage that it is not seriously 
affected by time variant coefficients. Certaine [6] has an approach similar 
to that of the present paper, but leaves much to be desired in effecting practical 
solutions. The methods given by Vasileva and Volosov [7-11] are promising 
but have received little use: They treat systems of differential equations which 

may be nonlinear and which include equations having small parameters as 
coefficients of the higher derivatives. 

The appendix of Andrus [1] contains a method for removing several 
large poles from a system of linear differential equations when it is known 
that terms in the solution, which correspond to these large poles, are insignifi- 
cant. 


Another approach which can effect some time savings involves partitioning 
the differential equations into two or more subsystems and integrating each 
subsystem with a different integration step-size depending upon the response 
time of the subsystem. 


THE NUMERICAL EFFECT OF SMALLTIME CONSTANTS 


Solutions to systems of linear differential equations with constant 
coefficients contain terms of the form ae^ where a and A are real or complex 
numbers. (The constant A is often called a pole or eigenvalue of the system, 
and the reciprocal of | A| is known as a time constan t when A is a negative real 
number.) Many numerical methods of integration are derived under the 
assumption that the solution can be expressed over an interval from T to T + At 
as a Maclaurin series 

b 0 + b t (At) + b 2 (At) 2 + ... 


truncated after several terms. However, if I A | is large, many terms of the 
series may be required to represent ae^ accurately. The Maclaurin expansion 
of e^ is 


A( t + At) At 
e = e 


1 + 


AAt 

i.» 


(AAt)' 

21 


Since | A| • At cannot greatly exceed unity in order for the series to be truncated 

after several terms and still approximate ; the maximum acceptable 

integration step-size is approximately equal to 1/ I A | , where I A I is the 

max max 

magnitude of the pole of largest absolute value. Even if the numerical method 
of integration is not based directly upon a Maclaurin expansion, the large 

magnitudes of higher derivatives of ae^ may still cause difficulty when |A| is 
large. 


Often the coefficient a is so small in magnitude that the term ae^* gives 
rise only to negligibly small variations in the solution (see Appendix A, Example 
2). However, if | A I is large, the higher derivatives of ae^ can still be so 
large in magnitude that one is forced to integrate with a very small integration 
step-size. 


DEFINING THE RESPONSE OF A TRANSFER FUNCTION BY MEANS 
OF A SYSTEM OF FIRST-ORDER LINEAR DIFFERENTIAL EQUATIONS 


Suppose that the transfer function L (E ou ^)/L(Ei n ) has been separated 
into partial fractions: 
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= C + 


( 2 ) 


L(E cut> 

L( V 


n 


Z 

i=l 


n. 


y.. 

__±L 


. . n.-j+l 

]=1 (s-A.) 1 


where c, yy, and Aj are constants and yy and A^ may be complex. The 

quantities Aj(i=l, 2, n) are the poles of the transfer function, and Aj is 

a pole of multiplicity nj. It will be shown that the following system of differential 
equations has the same transfer function: 


n. 

n i 



e = ce. + y. y. h.. 

out in f-' . j i] 

1=1 j=l 


(3a) 

h., - A.h = y E. 
il l il il in 


(3b) 

h. . - A. h. . = h. . , + (y. . - y. . . ) E. 
i] i ij 1,3-1 i] i,J-l m 

(j=2,3,. . . , n.) 

(3c) 

h..(0) =0 (j = 1,2, .... n.) 

i o 


(3d) 


where i = i , 2, . . . , n. 

On taking the Laplace transform of the members of (3b) and (3c) , one 
obtains : 


Y. 

L( h H )=— ^ L(E ) 
il s-A. m 


L(h ) = 


1 

L(h. 

1 9 

s-A. 

l 

, n.. 

Then 


y - v 

lj ~ L(E. ) 

s-A. in 

l 
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S L <v = 

j=i J 


r ii 


s-A. 


i J 


% 


Vg - y ii , T is ■ 


il s -V s <s -v ! 


r lk T i,k-1 

+ . . - + 2 — 


S-A. 


1 -I 


r ii 


( s -A. ) ‘ 


y i2 - y il 

s -\ 


i3j!i2l + J _hi 

^ J [<**i 

|_(s-A.) n i 


y i2 ~ y il 

/ -v \ n-i 

(s-A.) i 


y i2 ' Y il 


(s-A.) 


k-i 


y. -y. 
i n. i, n. . 
l l-l 

s-A. 

l 


L( E. ). 
m 


Simplification yields 


n. 

l 


r i'll 


Z Mb ) = 

j+i J 


i2 


/ ,, n - 

.(s-A.) i 


, . n.-l 

(s-A.) i 


7 in. “i 
1 


S - A i- 


L(E. ). 
in 


(4) 


The transfer function (2) may be easily derived from (3a) and (4), 


Since equations (3) have the transfer function (lc) , the solution to (3) 
is the solution defined by the transfer function. 

If the transfer function ( lc) has only real coefficients, then corresponding 
to the equations (3b) and (3c) , where Aj is any complex (non -real) pole, there 
are the redundant equations 

h. - A. h. = y. . • E. 
ll i ll il in 


h. . - A. h. . = h. . + (y . . - y. . , ) E. 

ij i ij i,J-i i] 1.3-1 m 


(j=2,3, . . . , n ) 


also present in ( 3b) and ( 3c) . Here the symbol x stands for the complex con- 
jugate of x. 

A more important observation is the following: If | A^l is large and A^ 

is negative ( in a stable system there can be no positive poles) , it is often true 

that the terms 11, corresponding to A., are negligibly small. (The relative 
xj i 

magnitudes of the constants y. (i= 1, 2, . . . , n; j = 1, 2, . . . , n.) may be the 
6 




deciding factors in such a case. However, the nature of E. is also important. ) 

If h is negligibly small for j = i, 2, . . . , n., then the corresponding differ- 
kj 1 

ential equations of (3) may be discarded and the remaining equations integrated 

numerically using a larger integration step-size. However, if the contributions 

of the terms h^ are significant, the numerical problem cannot be completely 

eliminated, because in such a case a small step-size is intrinsically necessary 
to determine the rapid responses of the solution E 


out’ 


But even in the latter 


case it is often possible to increase the step-size significantly by assuming some 

form from E over short time intervals and integrating equations (3) analyti- 
m 

cally rather than numerically. The remainder of this paper is concerned with 
the analytic solution to (3) and the manner in which it may be tied in with the 
numerical integration of the nonlinear equations of ( 1 ) . If the solution derived 
in this paper is employed, one need not concerned himself with eliminating 
negligible equations; the analytic solution essentially enables one to hop over 
negligibly small variations in E 


In passing, we remark that when| A.. | is large, the equation h.^-X.h.^y.^ E. q 
may be approximated in some cases by the equation h^= - (y.^/X.) E.^. This 
and similar approximations are discussed by Cohen [12]. 


REPRESENTATION OF E in AS A LINEAR FUNCTION OVER 
SHORT INTERVALS OF TIME 


Some form for E in must be assumed before equations (3) can be 
integrated. The function Ej n will be approximated over the time interval from 
T to T + At by means of the line passing through the two points 

[T, E in (T)], [T+ At, E in (T+ At)]. 

The corresponding expression for E^ n is 

E in (t) = E in (T) + [E in (T+ At) -E in (T)]. (5) 


Of course it would be possible to represent E in (t) in other ways. The linear 
expression has been chosen in order to keep the integration formulas to be 
derived as simple as possible. 
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ANALYTICAL INTEGRATION OF THE LINEAR EQUATIONS 


Assume that equations ( 3) have been integrated up to time T and that 
values of Ei n (T) and E^ n ( T^At) have been given. Approximating E in (t) by 
means of (5) over the time interval from T to T + At, we will derive an ex- 
pression for Eout at T + At in terms of Ej n (T) , Ej^CI^At) , and hjj(T) where 
i = 1 , 2, . . . , n and j = 1, 2, . . . , nj. Later this expression will be tied in 
with the numerical integration of equations (la) and (ib) . 


Specifically , we will integrate the equations 

r . 

h - A. h = y E. 
ll l li li in 


h.. - A. h.. = h. + (y.. - y. . . ) E. (j = 2,3, . . . , n ) 
i] l 13 1,3-1 ' 13 i, 3 -i m 1 

from time T to T + At in order to obtain hjj(Ti-At) in terms of Ej n ( T) , 
E in ( T+At) , h a ( T) , h i2 ( T), . . . , h in .( T) . 

For simplicity a new time variable r = t - T is introduced. Let 


( 6 ) 


E. 

in 


[r] = E 


m 


(T + t) 


and similarly for the functions E ou t and hjj. In other words, the variable r 
will have its origin at t = T. Then 


E. 

in 


[r] = E. 


in 



(E. [At] - E. [0]} . 
1 in in J 1 


(7) 


Letting t be the independent variable and taking the Laplace transforms of (6) 
and (7) , one arrives at the following equations: 




8L( V " h ii |0) '\ L< V =T il L(E in ) 

sL < h ,j) -y°i = L < h i,j-i> + y - r i, i - I > L( v 


4 . E. [At] - E. [ 0 ] 

L(E. ) = - E. [0] + — m m 

in s in s 


At 


(3 = 2,3, ... , n.) 
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Let 


r ° = E in f °] 


(8) 


r l = 


E. [At] - E. [0] 
in 1 in 


At 


After eliminating L(E^ n ) and employing (8) , one obtains 
i 


UhJ = 


ii s-A. 


{ h iit°J + Y il + 


s s 


L(h. .) =^— <L(h. . .) + h. .[ 0] + (y.. - y. . .) 

i] s-A. ^ 1,3-1 i] 1 J 13 '1,3-1 


r r 

+ _L 

s s 2 


( j = 2, 3, . . . , n.) 


Successive substitutions yields 


h ik [ °i , r . (y ik ' r i,k-i ) , r i (y ik -yk-i* 

■+■ -f \ 


^(hij) \ j-k+l j-k+1 

k=i l ( s-A.) J (s-A.) J s 


, .. 3-k+i 2 

(s-A.) J s 
1 


for j = l, 2, n. and y. = 0 . 

1 *0 


Assuming Aj + 0 and separating into partial fractions, one obtains 


L(h 


w> ■ I 


h ik'°l 


k=l | (s-X.) 


Pw +r,<T ik' V i.k-l ) 


r j-k+i 

E 


0-1 (-X.) j - k -“ +2 (s-X.)“ <-X.) J - k+1 s 

1 1 ' r 


+ r > (1 'ik" Y l,k-l ) l 


j-k+1 . 


E 


i-k-o;+ 2 


j-k+1 


<«=i ( -x 1 , 1 - k - a+3 (s-x 1 )“ (V- k * 2 . (-\y s 


i-k+l. 2 



The inverse Laplace transform is 


j 1 h [0] • r j " k e X f 


1J 


k=l 




j-k+i Q!-i A.T 

- 7 - T e 1 


0=1 (-x.,J- k -“ + 2 («-i ); ( -x l} >- k+1 


+I ^ y ik\k-1 ) 


3 y k+ * ( j-k-Q'+2)T Q ' *e\ T j-k+1 + t_ 




a=l (-X.) 3 “ k ‘ Q!+3 (Q ! -i)! (-A.) j_k+2 (-A.) j_k+1 . 


for A. * 0. 

l 


Evaluating h_ at t = At and rearranging, one obtains 


V Atl = £ 


, am / 3-1 A. At 
(At) e l 


j 

+ r ‘ ' Jl, 


J 

+ ri ^ 


(/3-1)! 


r js -i 

/ A ^-At 

- 7 

(At) e i 

u 

a= 0 

( -X.f~ a a 1 


1 

?-i /0 . , . .. 0! A. At 

V (0 -a) (At) e i 

Li 

2=0 (-A ) 

p-a+ 1 

^ O'! 


i + AL- 
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where 0 = j-k+1 and A^ + 0. Therefore, 


j 

h. .rAtl = Y, h. . , . r01 • a. + r 0 b.. + r< b.. 

ij L q =1 i» j-0+1 1 J 10 0 13 O 1 ijl 


0 ) 


where 


r a 


. . , /3 -1 A. At 
(At 3 K e 1 


i/3 0-1)! 

j 


b «0 ■ } 


A<At 




,-XiAt _ y 1 (-X,. 


At) 


Q! 


ce = 0 o'! 


ijl 


(T i, j-/3 + l” y i,j-0* 


A, At 


(-A.) 


0+1 


-0 


-A-iAt 




Qi =0 


1 (-X.At) 


Q! 


O'! 



r O' i 



-XiAt ^V 2 

l 

(-X.At) 

«■* 



a =0 J 

> 


( 10 a) 


2 

where X. =£ 0 and the summation I s nugatory when 0 = 1 . 

1 G! =0 

When |X At| n ^ is much less than unity, in order to avoid subtraction of nearly 
equal numbers during the computation of b.^ and b.^, one should employ the 


following power series expansions: 
■ x ‘ at - 1 


-i (-X.At)" « (-X.At)" 

Li a \ 


a=0 


a] 


<2 = 0 


li 



-X,At fc 1 <- X i At >“] 

+ (-A.At) 

-AjAt ^ 2 

(-A.At)" ' 

a! 

e - L 

oil 

a= 0 


a = 0 

- 


oo (-X At) a 

= L (< x - P ) — 

a = /3 + 1 


Equation (9) also holds for the case in which = 0. However, for this 
case the coefficients (10a) must be computed as follows: 


a. 


_ lAt* 


0-1 


10 (0-1)! 
j 

I 


J 

b ij° = Jt ± (Y i,j-0+l " Y i,j-0 ) 


(At) 

0! 


(10b) 


1 


b ijl ^ (Y i,j-0 + l Y i.j-0 ) 

0 -1 


. 1M 


0+1 


(0+1)! 


Substitution of expressions (8) into equation (9) yields 

j 


x r Atl = T, h. . , , [ 0] • a. + E. [ 0] • c + E. [ At] • c (11) 

ij 1 J j± i,j-/3+l l 10 m 1 1 ij0 in 1 xjl 


where 


b b. 

, ill „ _ Hi 

C nn ll.-jA • • » C - 


"ijO ijO At ’ ijl 
According to (3a) 


At 


n. 


n x 

E out l AtI = C E m ' At ' + £ 2 h Ji ( At > • 


( 12 ) 


i=l j=l 


12 



Equations (11) and (12) may be rewritten in terms of T as follows: 


n. 

n 1 

[ E out (T + At) = c E in (T + At) + S Z h (T + At) 

i=l j=l J 




j 


( 13 ) 


h.. (T + At) = Tj h. . (T) ' a + E. (T) * c..„ + E. (T+ At) • c... 

ij tL. i,j-/ 3 + 1' ' i (3 in ' ' ijO in' ' ijl 

P 1 


Equations (13) are the desired integration formulas. As long as At and the co- 
efficients of the transfer function (2) are constant, the coefficients c, a. , 

c.._, and c.._ will be constant. 
ijO ijl 


One can compute E , and 
out 

E . 
out 

from 

n 

n. 

l 


E (t) = cE. (t) + Yj 

out ' in ' 7 H. 

i=l 

z 

j=l 

h.. 

ij 

n 

n. 

i 


E (t) = cE. (t) + Z 

out m ' . ' 

i=l 

Z 

j=i 

h.. 

ij 


where h„ and h_ may be computed from equations (3) . 

In most problems n t = n 2 = . . . n = 1. Then (13) reduces to 


E out (T+ At) = cE. n (T + At) + 2 h n (T+At) 

i=l 


(T+ At) = h il (T) • a„+ (T) • c„ A + E^ (T + At) • c, 


where 


= b. 


b... 

ill 


il in 


ilO in 


ill 


c,.. = 


ill 


ilO ilO At ’ ill At 


13 



and where 


y il 


AjAt , 

a. = e 1 , b... = 

ll llO A. 

1 


(e X i at - l) 


b ill“ \. 2 


il . A^ At j , , 
(e - 1 - A^At) 


when A. # 0 and 

l 


a il 


b ilO ’ y il At ' b ill =y ii 


(At) 2 
2 


when A . = 0. 

i 


TIE-IN OF THE ANALYTICAL SOLUTION TO THE 
NUMERICAL INTEGRATION OF THE NONLINEAR EQUATIONS 


Formulas (13) may be tied in to the numerical integration of the non- 
linear equations 


I y k = V*’ y *' y2> ' • ■ • y „ d ’ E in- E out» (k * 2 V 

E in ’ * y >- y2 V ' V E „ut> 

d 


in a rather simple manner. Essentially E ^ may be eliminated from these 

equations by substituting the righthand member of the first of equations (13) for 
E Assuming the equations have been integrated up to time T, we may write 


n 


n. 


y k = 0 k Yl> y2 ’ ■ 


E. , cE. + Yj y, h..) 
in in i=l j=l « 


for t > T, where At = t - T and 


14 



h.. 

ij 



j-^+i 


(T) 


i/3 


E. (T) 
in 


c + 
ijO 


E. 


in 


c 


ijl’ 


and similarly for the equation E.^ = <p. Then the differential equations may be 
integrated numerically in the ordinary manner. 

If E is continuous , the approximate solution obtained by the method 
proposed in this paper will approach the exact solution as At approaches zero. 


THE CASE OF LINEAR EQUATIONS HAVING TIME 
DEPENDENT COEFFICIENTS 


If the coefficients of p(s) and g(s) in the transfer function (lc) vary with 
time, then c, A. and y will also be time dependent. This problem may be 

treated by the methods proposed, provided At is chosen sufficiently small so 

that c, A. and y. may be assumed to be constant over each At interval, 
i i] 


DERIVATION OF THE ANALYTICAL SOLUTION DIRECTLY 
FROM A SYSTEM OF LINEAR D IFFERENTIAL EQUATIONS 


Suppose that the relationship between and E ^ is expressed by means 
of the equations 


q 


Aq + E. b 
— m — 


"out 


T 

w E. + u q 
in — — 


( 14 ) 


rather than by a transfer function. Here A symbolizes an m x m matrix of 
constant elements, the symbols q, b, and u represent m x 1 vectors of constants, 
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T 

u is the transpose of u, and w is a scalar constant. It will now be shown 
briefly how equations (3) can be derived from (14) employing classical methods 
for the solution of linear systems. 


It is proven in elementary matrix theory that there exists a nonsingular 
matrix S such that 

S" 1 AS = J 


where J has the form: 


f Ji 0 . . . 0 ^ 

0 J 2 . . . o 


J = 


\° °"- V J 

where the 0's represent rectangular arrays of zero elements and 

^/p. 00... 000 

1 p. 0 . . . 0 0 0 

l 


J. = 

l 


0 1 p. 


0 0 0 


0 00... 01 p , 

V n; 


The constants p 1# p 2 , .... p^, are eigenvalues of A and are not necessarily 
distinct. The matrix J is known as the Jordan canonical matrix. 


Since A = SJS 1 


S -1 q = J(S" 1 q) + E. n (S^b) 

) T l 

E , = w E. + (u S) (S q) 
out in — 
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rp rp 

Defining £= S -1 a> £ = S _1 k> and v = u S, the above equations may be expressed 
as follows: 


£ = J £ + E ln S 


E = w E. + v £ 
out m ““ 


(15) 


Decompose the vectors £, g, and v in the manner indicated below: 

T x T T, 

£ ~ (£1 * £2 » • • • » Eq' ) » 


T T T T 

& = (Si » £2 Stf ) » 


T T T T\ 

v = (Vi , Y2 ) , 


where £., g., and v. (i = 1, 2, . . . , n^ are n' x 1 vectors. Now equations (15) 
may be broken down as follows: 

f% = V V E m% < i=1 ' 2 


_ y T 

= w • E. + 2 j Zj Ei • 

i=l 


Letting 


£. = (P ir P i2 > •••> P inr ) > 


% ^ g il’ g i2’ ' * * ’ g in^ ’ 


y. = (v. r v. 2 , .... v. n .) , 
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one may obtain the equations 


^il ^i ^il ^in ^il 


p.. - p. p.. = p. . , + E. g.. (j = 2, 3, . . . , n.> 

13 *1*13 *1,3-1 m i) r 


for i = 1, 2, . . . , n' . 

Now let f = v_p_ and multiply the differential equation defining p_ 
through by v_. Thus 

''ll' f ii ' E in (v il s il ) 

J V Vij = f i,j-l + E in (v ij S ii> (1=1 - 2 3 n I> 

n. 

n i 

E = w E. + Yj T, f.. (16) 

^ out in ij 

1=1 j=l 

Assuming £(0) = 0, it is obvious that f.. (0) = 0. The analytical solution of (16) 
may now be obtained in the same manner as the solution to equations (3). 

The transfer function of (16) is 


L(E OUt ) 

L <V 


< Z 


n l v.. g.. 

= w + £ Z k=1 * * , 

i=i j=l (s - P t ) n i 1 


Provided p t , p 2 , . . . , p are distinct, one may identify f.., p., w, n , n. and 
3 n lj i i 

E v., g, with h , A.., c, n, n., and-y.., respectively. However, in the rare 
. lk °ik ij i i 13 

k=l J 

event that the p.’s are not distinct, the identification is somewhat more corn- 
'll 

plicated. (For example, if p. = p and p. + p. for j #i, k, then max (m , il '") 

1 K 1 J 1 K 

must be identified with n. and n^ is zero. ) 
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T 

For the case n. = i, the vectors s., and r are simply right and left 

1 1 1 T 

eigenvectors of A corresponding to the eigenvalue such that r .^_s = 1. 

Therefore, for the usual case in which n = 1 for i = 1, 2, . . . , n, the coefficients 

y may be determined from the eigenvectors of A; in fact, y would be equal 
11 T T ll 

to (u _Sj^) • (r b) provided p. + p^ for i # k. If b has all zero elements except 

one — say the first element bj of b — then 

y n = b i , (17) 


where the first element of the eigenvector _r.^ is arbitrarily chosen to be unity. 

Equation (17) also establishes a relationship between the eigenvectors of A and 
the numerators of the partial fractions of the transfer function (2). 


CONCLUSIONS 

Equations (13) provide a simple integration formula for determining the 
response E at time T + At, of the transfer function (2) from values of the 

forcing function E. at times T and T + At. The section entitled "Tie-in of the 
in 

Analytical Solution to the Numerical Integration of the Nonlinear Equations" shows 

how equation (13) may be easily tied in with the numerical integration of other 

differential equations relating E. to E , 

m out 

The coefficients in (13) depend upon At and the coefficients of the transfer 
function in the manner shown in equations (10a) and (10b). Furthermore, the 
section entitled "Derivation of the Analytical Solution Directly from a System of 
Linear Differential Equations" derives the relationship of the coefficients in (13) 
to the eigenvalues and eigenvectors of the coefficient matrix A of the system (14) 
of first-order linear differential equations. (Indirectly, this also gives the de- 
pendence of the coefficients of the transfer function upon the eigenvalues and 
eigenvectors of A. ) The coefficients of (13) must be recomputed if either At or 
the coefficients of (2), or the alternate equations (14) have changed significantly. 

The integration formula (13) very often enables one to integrate numerically 
a system of control equations such as (1) at a much larger integration step-size 
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than the largest which could be used by conventional methods, especially when 
the transfer function has poles of large magnitude. 

Use may also be made of the first-order differential equations (3) which 
define E very concisely in terms of E (without introducing troublesome 

high derivatives of a single variable). If no difficulty with poles of large mag- 
nitude is present, equations (3) provide a simple means for determining E . 

O Ui 

They may be integrated numerically, simultaneously with other differential 
equations, and they may be solved analytically when E is available as an 

explicit function of time. When a pole A^ of large magnitude is present, it is 

often possible to eliminate it by simply removing from (3) the equations corre- 
sponding to Aj. However, a careful analysis of the particular problem should be 

made before any equations are eliminated. For this reason one may find it 
most convenient to utilize the integration formula (13) , which is likely to be a 
faster method even when poles of large magnitude are not present. 
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APPENDIX A 


EXAMPLES 


Example 1 

A problem in which the forcing function is given explicitly by 


E. = sincet 
m 


will now be examined. In this case equation (3 b) is 


il 


A.h.. = y.. sinat. 
1 xl il 


Assuming h. ( 0 ) = 0 and integrating, one obtains 


h (t) = — o o (ol e 1 - cxcos at - A. sinat) . 

il v ' + a i 


If A. < 0, 



Zii 

f 

a 

\ 

h il 1 s 

A. 

1 

l 2 

X, 

l 

V 


(18) 


In particular, consider the problem in which a = 1 and the transfer 
function is 


L(E )/L(E. ) 
out in 


4s 3 + 233s 2 + 9 98s + 5440 

2s 4 + 224s 3 + 2444s 2 + 440s + 4000 


(5/4) + 

s + l + \TT 


-(5/4) 'sTT 
s + 1 - si -1 


1 1 
+ , 
s + 10 s+100 
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Therefore 


X t = -1 -nTT , X 2 = -1 +nTT , \ 3 = -10 , X 4 


Tn = (5/4) \TT , y 21 = -(5/4) nTT, y 31 = 1 

It follows that 


hn (t) + h 21 (t) 


i [«'* 


Y41 


(2 cos t + sin t) - 2 cos t + sin 


in t j 


^31 (*) “ 


101 


(e - cos t + 10 sin t) 


1X41 (t) = 10001 ( e ~ 100t " cost+ 100 sint). 


-100 

1 . 


In Figure A-l, h n (t) + 

h 2i (t), h 31 (t) , and E ^ (t) are plotted as 

functions of t. Inequality (18) reveals that |h 41 | so. 0102. Therefore, the 
variable h 41 (t) is too small in magnitude to be plotted on the same scale, and 
yet the largeness of |X 4 | would make it necessary to use a numerical integration 
step-size roughly equal to 0. 01. A step-size of 0. 1 could suffice if the differential 
equations defining h 41 were eliminated. 


Example 2 


The electrical network shown in Figure A-2 is a filter used in an actual 
control system of a space vehicle. It has been represented by a system of nine 
differential equations having the for m indicated in (14). The response E ^ was 

determined by means of the integration formula (13) where E was defined by 

E. = sin (27rt). 
in 

Linear interpolation was used to approximate E over each At time interval. 

Time increments of At = 0. 2 sec, At = 0. 06 sec, and At = 0. 02 sec were 
employed. 

For comparison the integration was also performed by the fourth order 
Runge-Kutta method of numerical integration using At = 0. 006 sec. The latter 
technique diverged when an increment of At = 0. 007 sec was employed. The 
results of the comparison are shown in Figure A-3. 
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APPENDIX B 


MATHEMATICAL SUBROUTINES 

Subroutine for Decomposition of a Rational 
Fraction into Partial Fractions 


Input Data. 


A. (i= 1, 2,..., N r ) 
n. (i = 1, 2,. . . , N r ) 

N c 

Re (A ) , Im (A^) (i=l, 2, 


M 

a k (k= 0, 1, . . . ,M) 


number of real poles 
distinct real poles 

multiplicities of real poles (usually n = 1) 
one-half the number of complex poles 
, N ) real and imaginary parts of complex poles 

v 

degree of numerator of rational fraction 
coefficients of numerator 


Note : A pole is defined to be a root of the denominator of the rational fraction. 

All complex (non-real) poles are assumed to be distinct. The conjugate 
of each complex pole given in the input data is also a pole. 

The rational fraction must have a numerator of smaller degree than the 
denominator. The coefficient of the term of highest degree in the de- 
nominator is assumed to be unity. 

Definition. We define the polynomials 
M M-l 

P(s) = a o s + a i s + • • • + a M _i S+ a M 
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(i=i,2,.. ,N r ) 


^(s) = 


' N 



n R (s-x) n k 

• 

n C (s 2 + b s + c ) 

= 1 


_ k = 1 


(k*i) 


where 


b k = -2Be ( A k ), o k =[Re<\)] 2 + [ Im <v] 
Computation . Compute 


P(*j) 

7 il ~ h. (X.) 

1 V 1 


y.. = 


r d ^- 1 


ij (]-!)•' 


P(s) 

dJ " 1 IV S) 


(j = 2, 3, . . . , n.) 


s = X 


for i = 1, 2 N r . 


Compute the complex numbers 

p(X.) (i=l,2, . . ,N C ) 



• 

n c < A i + Vi + V 

_ k = 1 


.(k*i) 


where 


A. = Re(A.) + \T- Tlm(A.) (i = 1, 2, N^) are complex. 

Ill L 

Compute 

Im(Q.) 

r il “ Im(AA ’ 
for i = 1, 2, . . . , N c> 


r i2 = _r il ' Re( V + Re(C V 
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Output Data . 


Y ij (i = 2 > • • • . n r 5 j = 1, 2, . 

r u- r 2i < 1 = 1 - 2 V 



coefficients of 
numerators of 
partial fractions 


Analysis. The procedure described above decomposes the rational fraction 


M , M-l 

s + a t s + 


+ a M-l S + a M 

n r 

n (s - a ) k 

k = 1 k 


'N c 

n (s 2 + b s + c ) 

k = 1 k k 


into partial fractions: 
N_ n 


R 

Z 

k=l 


z k 

j=l 


y 




(s 


V 


n k - j + 1 


N c 

- z 

k=l 


r kl S+r k2 
s 2 + b, s + c. 


( 1 ) 


It is assumed that the rqots of s 2 + b^s + are the given complex numbers 

and A^ conjugate. It is also assumed that the numerator of (1) has lower degree 

than the denominator, that the real roots A 2 , . . . , A are distinct, and 

W R 


that there are no repeated complex roots. The method of decomposition is 
standard. 


Subroutine for Computing Coefficients Defining the 
Response of a Transfer Function 


Input Data . 

N number of real poles 

X\ 

A. (i = 1, 2, . . . , N ) distinct real poles 

l H 
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multiplicities of distinct real poles 


n. (i = 1, 2, . . . , n.) 


y. . (i=l, 2, . . . , N • j=l, 2, . . , n.) numerators of partial fractions 

1 j -H 1 X 7 


N, 


corresponding to real poles 
one-half the number of complex poles 


Be (A ) ,Im(A ) (i=l, 2, . . , N c > real and imaginary parts of complex 

11 ^ poles (It is assumed that there are no 

repeated complex poles. ) 


hi' r i2 <‘ =1 ' 2 N C> 


coefficients of numerators of partial 
fractions corresponding to complex poles 


At 


time interval 


Computations . Compute 


a.. = 
ij 


j = l >2,. • ,n 


. 4 A. At 
( At' 1-1 e 1 

(j " !)•' 


h = V Y Ail-jS+l y . h J zM 

iiO ^ R 

J (8=1 (~X.) P 


A. At (8-1 (-A.At) 
l-e 1 l 1 


a 


a=0 


a l 


iii ^ <v +1 ' 


A.At (3 -1 (-A.At) 

l-e 1 Yj 1 


a 


a=0 


a: 


+ (-A.At) 


A.At d-2 (-A.At) 

l-e 1 E ^ 

LJ ctl 


(X 


a=0 


fc 2 


assuming A =£ 0, where y. _ is defined to be zero. The sum Y. is zero when 

10 a = 0 

(9* »• 
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n. 

If IAj • At | < 0. 01, then the term in brackets in the expression for b should 

be computed by means of the series 


A. A °° (-A At) 

e‘ Z 


a 


a = j3 


a: 


and the term in braces is the expression for b ^ should be computed by means 
of the series 


A. At oo ( -A. At) 

>* Z 


a 


a = p + 1 


a! 


■ (a-p) 


When A^ = 0, compute 


a.. = ^ 


3-1 


ij (3-1)! 
j 


b iJ0 = J =1 <r i,J-0-n"W 
j 

^ b ijl Jii (y i,j -J8-1 " y i,j-/3 ) 


m n 


p- 


(At) 


JS+1 


W + 1 )! 


where, again , y = 0. 
Also compute 


c = b 
ijO ijO At 


, = J*L 

ijl At 


for i = 1, 2, . . . , N and j = 1, 2 n.. 

it 1 
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Compute 


Re (A.) • At 


Re (A^) = e 


cos Im ( A.) • At (Note: = e 


A. At 
1 


Im(A.) 


Re (A.) • At 
e 1 sin 


in Jim (A. • At) j 


Re(r') 


r i/ 2 


Im(I\ ) 


B 


iO 


T. 0 + T. . • Re (A.) 
i2 ll x r 

2 • Im(A.) 

r.' 

±- ‘V 1 ’ 


} <i=l,2,..,N c ) 


B 


il 


r. 

= ^2- (A. - 1 + A. • At) 
i 


c i0 


= B 


'il 


iO At 


c 

il At 


where B.„, B.„, rT , A., A., C.„, C.., are complex numbers. 
iO il i l l iO il 

If |A^At| <0.01, compute 


i. - 1 = - A. Yj 

1 l j 

a-1 


- (-A.At) 


a 


oil 


i. - i + A.At = - A. T] 

l l l ^ 


S (-A.At) 


O' 


a =2 


oil 


by series as indicated. 
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Output Data . 


coefficients defining 
the response of a 
transfer function 


a.. 

i] 


ijO } 


■N 

i= 1> 2, 

N 

> 

j = 1, 2, .. 

n. 


C iji 


Re (A.), Im(A.) 
Re(C. 0 ) 5 Im(C i0 ) l 
Re^), ImfC.^ 


i - 1, 2, 







Analysis. The derivations of the coefficients are contained in the main part of 
the paper. The coefficients are for use in a subroutine for determining the re- 
sponse of a transfer function of the form 


L(E .) 
out 

L(E. ) 

' in' 

Subroutine for Determining the Response 
of a Transfer Function F(s) 


N. 


R 


n, 


k 


= c + 


_kL 


2 2 „ 

k=l (s-A ) n k 3 


N C 

- Z 

k=l 


r kl S+r k2 
s 2 + b, s + c. 


Input Data . 

N number of distinct real poles 

X\ 

one half the number of complex poles 
n. (i = 1, 2 N r ) multiplicities of real poles (usually n. = 1) 
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c (c is usually zero) 


V C ijO’ 

O 

w* 

t— 

II 

2,.. . ,N r ; j=l,2,.. 

•,n.) 



Re (A.), 

Im(A.) 

(i=l, 2, . . 

N c> 


) 

coefficients 

R e (C i0 ), 

^(ao) 

(1=1,2,... , N c ) 




RefC^), 

Im/C.j) 

(i =1 , 2, • . 

N c> 





E. (T) 
in' ' 

h Aj (T) (i 
Re[H.(T)] 


E. (t) 
in' ' 


Computations . 


=1,2 N r ; 

, Im[H.(T)] 


j = 1, 2 n.) 

(1=1,2 N c ) 


initializing 

conditions 


current value of forcing function 


n r n i 


N, 


where 


E (t) = c- E (t) + Yj Y h (t) + 2 Y Re[ H (t)] 
OUt in i=l j = l 1J i=l 


J 

h..(t) = Y h- • „ , 4 ( t ) -a. + E. (T) • c + E. (t) • c.. 
lj ’ “ l, J-/S+1 i/3 in' ' ijO in' ij 

P 


H,(t) = H l( T) • A,+ E |n (T) ■ C. 0 + E. n (t) • C. ( . 


ijl 


i m 


il 


Here H^/t), (T) , A^, C^, and C.^ are all complex. 


Note on H andling of Input Data . Since the subroutine may be used in one program 
to compute the responses of many linear components of a control system, and 
since the coefficients may change from time to time for any given component, it 
is probably advisable to place each possible set of input data into an area of 
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common storage. Then, at the time the subroutine is called, the desired set of 
data can be indicated in some manner. 

Even if the form of F(s) is left unchanged, more than one set of coefficients for 
the component may be required, because the coefficients are dependent upon the 
time interval At = T - t. And even if the basic integration step-size is constant, 
some methods of numerical integration will require evaluation of E j.(t) at 

fractional steps. (For fractional steps the computation of h..(t) and H (t) is 

lj k 

incidental, and it is not required that these values be stored for future use. ) 
Output Data . 

E (t) Response at time t 

h (t) (i = 1 , 2 , . . . , N ; j=i, 2, . . . , n ) 

IJ Hi 

Re[H.(t)], ImfH.(t)] (i=l, 2, . . ,N C ) 


new initializing 
conditions 


Analysis . Consider a transfer function 


L < E out> 

L <V 


c + f(s) 


( 1 ) 


where c is a constant and f(s) is a rational fraction in which the denominator 

has larger degree than the degree of the numerator. It is assumed that, before 

the subroutine has been called, the time histories of E. and E , have been 

m out 

determined up to time T and that the value of E. (t) for some time t > T has 

m 

been specified. Making the approximation that E is linear from time T to 

time t, one may determine approximately, using this subroutine, the corresponding 
value E j.(t) of the response. The derivation of the method employed is con- 
tained in the paper, which also explains how this subroutine may be tied in very 
simply with the numerical integration of differential equations other than those 
represented by (1) . 

The initial conditions at time T are expressed by means of the input variables 

E. (T) , h.. (T) , and H (T) . The behavior of f (s) is expressed by means of 
m ij k 
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the coefficients given in the input data. The assumption has been made that none 
of the complex (non-real) poles of f (s) are repeated. 

This subroutine can be used to determine the solution defined by any transfer 
function (1) with non-repeated complex poles. However, the primary purpose 
of the subroutine is to avoid the prohibitively small integration step-sizes re- 
quired by ordinary methods of numerical integration in the presence of poles of 
large magnitude. 
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